Nonlinear Registration

Introduction

This exercise involved implementing a nonlinear registration technique using the Gauss-Newton algorithm to minimize the sum of squared intensity differences. The target was to align a 2D slice from MRI brain scans of two different subjects. To treat world and voxel coordinates as equivalent, the images had already undergone preprocessing, which included alignment and affine resampling in 3D. In cases where tissue deformation needs to be addressed, linear functions used for fine registration are extended to nonlinear ones. Typically, any global disparities in scale and orientation have been corrected through previous affine registration, leaving only local deformations, denoted as $\delta$, to be adjusted:

$$ y_d(\mathbf{x}, \mathbf{w_d}) = x_d + \delta_d(\mathbf{x}, \mathbf{w_d}), $$ where $$ \delta_d(\mathbf{x}, \mathbf{w_d}) = \sum_{m=0}^{M-1} w_{d,m} \phi_m(\mathbf{x}) $$

and where $ \phi_m(\mathbf{x}) $ represents a set of $M$ basis functions, with $w_d = (w_{d0}, w_{d1}, \ldots, w_{dM-1})^T$ as the weights.

Nonlinear registration is a cruical tool in medical image processing for correcting distortions and geometric discrepancies between images. This process involves using basis functions to model the nonlinear transformations and adjusting the parameters to reduce the intensity difference between the fixed and moving images.

MRI images do not have quantitative physical intensity units, then preprocessing is needed to make the intensity values comparable. This can be achieved by minimizing the sum of squared differences between the intensities of corresponding points in the two images:

$$ E(\mathbf{w}) = \sum_{n=1}^N \left[ \mathcal{F}(\mathbf{x}_n) - \mathcal{M}(\mathbf{y}(\mathbf{x}_n, \mathbf{w})) \right]^2. $$

where $x_n$ denotes the voxel coordinates in the fixed image, $\mathcal{F}(\mathbf{x}_n)$ represents the corresponding intensity values, and $\mathcal{M}(\mathbf{y}(\mathbf{x}_n, \mathbf{w}))$ refers to the interpolated intensities in the moving image at the transformed positions $\mathbf{y}(\mathbf{x}_n, \mathbf{w})$. If the transformed positions lie outside the image boundaries, a default intensity value is used.

With transformations that have many parameters, such as nonlinear deformations, numerical optimization methods are necessary. The Gauss-Newton algorithm is well-suited for this, as it linearizes $\mathcal{M}(\mathbf{y}(\mathbf{x}n, \mathbf{w}))$ with respect to the weights $w$. The derivative with respect to $w{dm}$ is given by:

$$ \frac{\partial M(y(x_n, w))}{\partial w_{dm}} = g_d^n \cdot \phi_m(x_n). $$

where $g_d^n$ is the partial derivative of the moving image along dimension $d$. The optimization process uses a Taylor series approximation to iteratively update the parameters. If the update step becomes too large, the Levenberg-Marquardt method is applied to maintain stability.

The main aim of this exercise is to correct remaining anatomical differences by minimizing intensity differences between corresponding points in the two images.

Input data and code hints

Import Python libraries:

import numpy as np
np.set_printoptions( suppress=True )
from matplotlib import pyplot as plt
plt.ion()
from scipy.ndimage import map_coordinates
import nibabel as nib

Read the two 3D scans you'll be working with in this exercise and select 2D slices:

fixed = nib.load( 'IXI092-HH-1436-T1.nii.gz' ).get_fdata()
moving = nib.load( 'IXI097-HH-1619-T1_resampled.nii.gz' ).get_fdata()

fixed, moving = fixed[ ..., 73 ], moving[ ..., 73 ] # get slices

Below are three functions which let you visualize the images and corresponding grids in the same way. You can use these functions through all tasks where needed.

def visualize_image( ax, im, **kwargs ):
    # Visualize an image the same way the grid and vector fields are drawn: to the right means
    # x1 (row number in image matrix) goes up, and to the top means x2 (column number in image 
    # matrix) goes up.
    ax.imshow( np.flipud( im.T ), **kwargs )
    ax.set_aspect( 'equal' )
    ax.axis( False )

def plot_grid( ax, X1, X2, **kwargs ):
    for i in range(X1.shape[0]):
        ax.plot(X1[i,:], X2[i,:], **kwargs)
    for i in range(X1.shape[1]):
        ax.plot(X1[:,i], X2[:,i], **kwargs)

def visualize_grid( ax, X1, X2 ):
    N1, N2 = X1.shape
    plot_grid( ax, X1[ ::10, ::10 ], X2[ ::10, ::10 ], color='b' )
    ax.set_aspect( 'equal' )
    ax.axis( False )

fixed = np.flipud( fixed )
moving = np.flipud( moving )

fig, ax = plt.subplots( 1, 2 )
visualize_image( ax[0], fixed, cmap='gray' )
visualize_image( ax[1], moving, cmap='gray' )
ax[0].set_title("Fixed Image")
ax[1].set_title("Moving Image")

X1, X2 = np.meshgrid( np.arange(fixed.shape[0]), np.arange(fixed.shape[1]), indexing='ij' )
fig, ax = plt.subplots()
visualize_grid( ax, X1, X2 )
fig.suptitle("Image Grid")

Text(0.5, 0.98, 'Image Grid')

In this exercise, you'll be asked to compute and display image histograms, and to indicate the location of a specific intensity value in the histogram plot. The code snippet below shows how this can be done:

data = fixed
value = np.median( data )

histogram, edges = np.histogram( data.ravel(), bins=100 )
binCenters = ( edges[ 1: ] + edges[ :-1 ] ) / 2
binWidth = edges[ 1 ] - edges[ 0 ]
fig, ax = plt.subplots()
ax.bar( binCenters, histogram, width=binWidth )
ax.set_ylim( 0, histogram[1:].max() )  # clip the entry for background
ax.set_xlim( 0, data.max() )
ax.plot( [ value ]*2, ax.get_ylim(), 'r' )
ax.grid()

Task 1: Normalize the intensity values of both images

Since MRI images are not calibrated, the pixel intensities in the fixed and the moving image are not necessarily comparable. As a first attempt to normalize the intensity values, divide each image by its maximum intensity value (to bring both images in the range $[0,1]$) and display both images, as well as their histograms, side by side. Are the intensity values now comparable between the two images?

Compute a robust version of the maximum intensity in each image by locating the intensity value in the histogram so that 99% of the image pixels have a lower intensity. For each image, indicate the location of the value you find in the histogram using a verticle line (see code snippet given in the introduction).

Now rescale the image intensities again, this time using the robust maximum value instead. (Intensities higher than the robust maximum should be clipped so that all intensities fall within the range $[0,1]$ after intensity normalization.) Display both normalized images, as well as their histograms, side by side again. Decide which normalization method makes the image intensities more comparable, explain why, and perform the remainder of the exercise using the normalized images of your choice.

import matplotlib.pyplot as plt

fixed_max = np.max(fixed)
moving_max = np.max(moving)

# Normalization
fixed_normalized = fixed/fixed_max
moving_normalized = moving/moving_max

print(fixed_max)
print(moving_max)

fig, ax = plt.subplots(2, 2, figsize=(10, 10))  # 2 rows, 2 columns

ax = ax.ravel()

visualize_image(ax[0], fixed_normalized, cmap='gray', vmin=0, vmax=1)
visualize_image(ax[1], moving_normalized, cmap='gray', vmin=0, vmax=1)
ax[0].set_title("Normalized Fixed Image")
ax[1].set_title("Normalized Moving Image")

# Histograms
ax[2].hist(fixed_normalized.ravel(), bins=256, color='gray')
ax[3].hist(moving_normalized.ravel(), bins=256, color='gray')
ax[2].set_title("Fixed Image Histogram")
ax[3].set_title("Moving Image Histogram")

plt.tight_layout()
plt.show()

795.0
1243.55224609375

# Compute the 99th percentile for both fixed and moving images
fixed_robust_max = np.percentile(fixed, 99)
moving_robust_max = np.percentile(moving, 99)

# Normalize and clip the images using the robust maximum
fixed_robust_normalized = np.clip(fixed / fixed_robust_max, 0, 1)
moving_robust_normalized = np.clip(moving / moving_robust_max, 0, 1)

# Plot histograms and mark the 99th percentile
def plot_histogram_with_robust_max(data, robust_max, title):
    histogram, edges = np.histogram(data.ravel(), bins=100)
    bin_centers = (edges[1:] + edges[:-1]) / 2
    bin_width = edges[1] - edges[0]
    
    fig, ax = plt.subplots(figsize=(6, 4))
    ax.bar(bin_centers, histogram, width=bin_width, color='gray')
    
    # Set the y-limit to avoid extreme outliers skewing the plot
    ax.set_ylim(0, histogram[1:].max())
    # Set the x-limit to the maximum intensity of the image
    ax.set_xlim(0, data.max())
    
    # Plot a vertical red line for the robust maximum (99th percentile)
    ax.plot([robust_max]*2, ax.get_ylim(), 'r', label="99th Percentile")
    
    ax.set_title(title)
    ax.grid()
    ax.legend()
    plt.show()

# Display images and histograms with robust max marked
plot_histogram_with_robust_max(fixed, fixed_robust_max, title="Fixed Image Histogram (99th Percentile)")
plot_histogram_with_robust_max(moving, moving_robust_max, title="Moving Image Histogram (99th Percentile)")

# Display the robust normalized images
fig, ax = plt.subplots(1, 2, figsize=(10, 5))
visualize_image(ax[0], fixed_robust_normalized, cmap='gray')
visualize_image(ax[1], moving_robust_normalized, cmap='gray')
ax[0].set_title("Robust Normalized Fixed Image")
ax[1].set_title("Robust Normalized Moving Image")
plt.tight_layout()
plt.show()

In the first Approach (Normalization by Maximum Intensity) the images are normalized by dividing each by its respective maximum intensity value. This scales both the fixed and moving images into the range [0,1] and while the images are now in the same intensity range, their pixel intensity distributions are not comparable due to differences in contrast. The histograms reveal that each image has a different distribution of pixel intensities, which means that direct comparison is still difficult.

In the second Approach (Robust Normalization by 99th Percentile) the intensity normalization is done using the 99th percentile value instead of the absolute maximum. We marked 99th percentile threshold just so we can ignore extreme outliers that could skew the results.Looking at the histograms we can see that this method provides a more robust intensity scaling by focusing on the majority of the pixel values (99%) rather than using the absolute maximum.

Task 2: Create basis functions for nonlinear deformation

In order to model the residual deformation in each direction $d=1,\ldots,D$, you are going to use a weighted set of basis functions: $$ \delta_d(\mathbf{x}, \mathbf{w}_d) = \sum_{m=0}^{M-1} w_{d,m}\phi_m (\mathbf{x}) . $$ Here $\phi_m (\mathbf{x})$ are $M$ basis functions, and $\mathbf{w}_d$ are the weights for direction $d$. For the purpose of this exercise, construct a set of $M=7 \times 7=49$ separable, 2D cubic B-spline basis functions that cover the image area of the fixed image. Include a visualization of the basis functions, plotted in a $7 \times 7$ grid, in your report.

For the purpose of the exercise, retain only the central part of $5 \times 5$ basis functions, removing all the B-splines that are centered on one of the image edges. This will avoid strangle-looking, large deformations of the image edges later on.

# cubic B-spline interpolation
def eval_cubic_BSpline(x):
    """
    Evaluates the uniform B-spline of cubic order "3" at the locations in vector "x"
    """

    x = np.asarray(x)  # Ensure x is a NumPy array
    y = np.zeros_like(x)

    # evaluate the uniform ubic B-spline at the locations in x
    abs_x = np.abs(x)
    
    # Different cases for cubic B-spline
    y[abs_x < 1] = (4 - 6 * abs_x[abs_x < 1]**2 + 3 * abs_x[abs_x < 1]**3) / 6.0
    y[(abs_x >= 1) & (abs_x < 2)] = (2 - abs_x[(abs_x >= 1) & (abs_x < 2)])**3 / 6.0
    y[abs_x >= 2] = 0

    return y

# Grid and spline setup
N_dim = 256  # Image dimension
M_dim = 7    # Number of basis functions in each direction

def generate_splines(N_dim, M_dim):
    # Knot spacing for the B-splines
    h1 = (N_dim - 1) / (M_dim - 1)

    # Create the B-spline basis functions for each dimension
    phi1 = np.zeros((N_dim, M_dim))
    phi2 = np.zeros((N_dim, M_dim))

    x = np.arange(N_dim)
    for i in range(M_dim):
        centered = i * h1
        shifted = (x - centered) / h1
        phi1[:, i] = eval_cubic_BSpline(shifted)
        phi2[:, i] = eval_cubic_BSpline(shifted)

    # Form the 2D B-spline basis functions using the Kronecker product
    return np.kron(phi2, phi1), phi1, phi2  # Full 7x7 grid of basis functions

Phi, phi1, phi2 = generate_splines(N_dim, M_dim)

# Now retain only the central 5x5 basis functions
Phi_central = np.kron(phi2[:, 1:-1], phi1[:, 1:-1])  # 5x5 grid of basis functions

# Visualize the central 25 basis functions in a 5x5 grid
fig, axs = plt.subplots(5, 5, figsize=(10, 10))

for i in range(25):
    ax = axs[i // 5, i % 5]
    basis_fn = Phi_central[:, i].reshape(N_dim, N_dim, order='F')
    
    ax.imshow(basis_fn, cmap='gray')
    ax.set_title(f'B-spline {i + 1}')
    ax.axis('off')

plt.tight_layout()
plt.show()

By visualizing the basis functions, we can see how each function will influence a specific region of the image. When we later combine them with the weights computed froBy visualizing the basis functions, we can see how each function will influence a specific region of the image. When we later combine them with the weights computed from the image, this visualization helps understand how local intensity values are smoothed or adjusted.

The result is a set of smooth 2D basis functions that can be used to model local medical image deformations, avoiding abrupt changes near the image edges.m the image, this visualization helps understand how local intensity values are smoothed or adjusted.

Task 3: Manipulate the image locally with basis function

In this task, you will use two weight vectors $\mathbf{w}_1$ and $\mathbf{w}_2$ to control the transformation model $$ y_d(\mathbf{x}, \mathbf{w}_d) = x_d + \delta_d(\mathbf{x}, \mathbf{w}_d), \quad d=1,2 , $$ and visualize the effect of changing these weights on the resulting image deformation. For each choice of weights, you should show (1) how a regular image grid defined in the space of the fixed image looks after deformation, and (2) the corresponding interpolated moving image, resampled to the image grid of the fixed image.

In a first experiment, set all elements of $\mathbf{w}_2$ to zero, and set all elements except one in $\mathbf{w}_1$ to zero. Try a few different values for the one non-zero entry, show each time the resulting deformated image grid and the interpolated moving image, and comment on what you see.

In a second experiment, apply a non-zero value to the corresponding (i.e., controlling the same 2D basis function) element in $\mathbf{w}_2$ as well. Again visualize the result and comment.

def delta_d(w_d, Phi):
    """
    Compute the deformation delta_d for each direction d.
    Arguments:
    - w_d: weights vector for direction d
    - Phi: 2D basis functions
    
    Returns:
    - The deformation for each point in the grid
    """
    return np.dot(Phi, w_d)

# Define the image grid and basis functions (assuming Phi_central from Task 2)
N_dim = 256  # Image dimension
n_bsplines = 25   # Number of central basis functions

# Create a meshgrid for the image grid
# X1, X2 = np.meshgrid(np.arange(N_dim), np.arange(N_dim), indexing='ij')
X1, X2 = np.meshgrid( np.arange(fixed.shape[0]), np.arange(fixed.shape[1]), indexing='ij' )


# Vectorize the grid for easier manipulation
grid_x = np.vstack([X1.ravel(), X2.ravel()]).T

# Experiment 1: Set w2 to zero and use a single non-zero value in w1
w1 = np.zeros(n_bsplines)
w2 = np.zeros(n_bsplines)

# Select one basis function to control the deformation
index = 12  # You can experiment with different indices
w1[index] = 50.0  # Modify this value for different deformation magnitudes

values = [25, 50, 75, 100]

# Compute the deformation for the different values of w1
for v in values:
    w1[index] = v
    delta_x1 = delta_d(w1, Phi_central)
    delta_x2 = delta_d(w2, Phi_central)
    
    # Apply the deformation to the grid
    deformed_grid_x1 = grid_x[:, 0] + delta_x1
    deformed_grid_x2 = grid_x[:, 1] + delta_x2
    
    # Reshape the deformed grid for visualization
    deformed_X1 = deformed_grid_x1.reshape(N_dim, N_dim)
    deformed_X2 = deformed_grid_x2.reshape(N_dim, N_dim)

    # Interpolate the moving image using the deformed grid
    def_interpolated_moving_img = map_coordinates(moving_robust_normalized, (deformed_grid_x1, deformed_grid_x2), order=3, mode='nearest')
    def_interpolated_moving_img = def_interpolated_moving_img.reshape(N_dim, N_dim)
    
    # Visualize the deformed grid
    fig, ax = plt.subplots(1, 2, figsize=(10, 5))
    visualize_grid(ax[0], deformed_X1, deformed_X2)
    ax[0].set_title(f"Deformed Grid 1 (w1={v})")
    
    # Visualize the interpolated moving image
    visualize_image(ax[1], def_interpolated_moving_img, cmap='gray')
    ax[1].set_title(f"Interpolated Moving Image 1 (w1={v})")
    
    plt.tight_layout()
    plt.show()

# Experiment 2: Apply a non-zero value to the corresponding element in w2
w2[index] = 50.0  # Control the same basis function in both directions

# Compute the deformation for the different values of w2
for v in values:
    w2[index] = v
    delta_x1 = delta_d(w1, Phi_central)
    delta_x2 = delta_d(w2, Phi_central)

    # Apply the deformation to the grid
    deformed_grid_x1 = grid_x[:, 0] + delta_x1
    deformed_grid_x2 = grid_x[:, 1] + delta_x2

    # Reshape the deformed grid for visualization
    deformed_X1 = deformed_grid_x1.reshape(N_dim, N_dim)
    deformed_X2 = deformed_grid_x2.reshape(N_dim, N_dim)

    # Interpolate the moving image using the deformed grid
    def_interpolated_moving_img = map_coordinates(moving_robust_normalized, (deformed_grid_x1, deformed_grid_x2), order=3, mode='nearest')
    def_interpolated_moving_img = def_interpolated_moving_img.reshape(N_dim, N_dim)

    # Visualize the deformed grid and moving image for Experiment 2
    fig, ax = plt.subplots(1, 2, figsize=(10, 5))

    # Visualize the deformed grid
    visualize_grid(ax[0], deformed_X1, deformed_X2)
    ax[0].set_title(f"Deformed Grid 2 w1 = 100 (w2={v})")

    # Visualize the interpolated moving image
    visualize_image(ax[1], def_interpolated_moving_img, cmap='gray')
    ax[1].set_title(f"Interpolated Moving Image 2 w1 = 100 (w2={v})")

    plt.tight_layout()
    plt.show()

In this task, the goal was to experiment with different weights to control the deformation model applied to an image grid. First, we set all elements of the weight vector W2 to zero and chose one non-zero element in W1 to see how changes in this single weight affect the image grid and the deformed moving image. We then applied the deformation to the grid, and the moving image was resampled according to the new grid structure. As the single weight in W1 increased, the grid lines in the corresponding region were pulled more strongly, causing visible localized deformations in both the grid and the resampled moving image.

In the second experiment, the same weight index in W2 was set to non-zero values, introducing a 2D deformation by controlling the same basis function in both directions. This added more complex, two-directional grid deformations, resulting in a greater shift and stretching of the grid and the image. The results showed that applying non-zero weights in both W1 and W2 allowed for more dynamic and flexible deformations compared to controlling only one direction, resulting in a broader range of image distortions and transformations.

Task 4: Implement Gauss-Newton optimization

Search for the weight vector $\mathbf{w} = (\mathbf{w}_1, \mathbf{w}_2)^T$ that minimizes the sum-of-squared-differences energy $$ E(\mathbf{w}) = \sum_{n=1}^N \left[ \mathcal{F}(\mathbf{x}_n) - \mathcal{M}(\mathbf{y}(\mathbf{x}_n, \mathbf{w})) \right]^2 , $$ using the Gauss-Newton algorithm described on page 24-25 in the book. Start from a zero weight vector $\mathbf{w}$ (corresponding to no deformation), and update $\mathbf{w}$ for 100 iterations. For the first and every 20-th subsequent iteration, show both the deformed image grid and the resampled moving image. Also plot how the energy changes across all the 100 iterations.

# Parameters
N_dim = 256  # Image dimension
n_bsplines = 25   # Number of central basis functions
delta = 0.001  # Small step for finite differences
iterations = 100  # Number of iterations

# Initialize weight vectors
w1 = np.zeros(n_bsplines)
w2 = np.zeros(n_bsplines)

# Initialize energy array to store energy at each iteration
energy_history = []

# Gauss-Newton optimization loop
for iteration in range(iterations):
    # Compute the mapped coordinates
    delta_x1 = delta_d(w1, Phi_central)
    delta_x2 = delta_d(w2, Phi_central)
    deformed_grid_x1 = grid_x[:, 0] + delta_x1
    deformed_grid_x2 = grid_x[:, 1] + delta_x2
    
    # Resample the moving image using the mapped coordinates
    moving_resampled = map_coordinates(moving, [deformed_grid_x1, deformed_grid_x2], order=1)

    # Compute the residual (difference between fixed and resampled moving image)
    residual = fixed.ravel() - moving_resampled
    
    # Compute the energy (sum of squared differences)
    energy = np.sum(residual**2)
    energy_history.append(energy)

    # Compute Jacobian matrix as described in the slides
    G1 = (map_coordinates(moving, [deformed_grid_x1 + delta, deformed_grid_x2], order=1) - moving_resampled) / delta
    G2 = (map_coordinates(moving, [deformed_grid_x1, deformed_grid_x2 + delta], order=1) - moving_resampled) / delta

    # Jacobian: J = [G1 * Phi_central, G2 * Phi_central]
    J = np.hstack([G1[:, np.newaxis] * Phi_central, G2[:, np.newaxis] * Phi_central])

    # Solve the linear regression problem: (J^T J)^{-1} J^T r
    JTJ = J.T @ J
    JTr = J.T @ residual
    update = np.linalg.solve(JTJ, JTr)
    
    # Update weights
    w1 -= update[:n_bsplines]
    w2 -= update[n_bsplines:]
    
    # Visualization for the first and every 20th iteration
    if iteration == 0 or (iteration + 1) % 20 == 0:
        fig, ax = plt.subplots(1, 2, figsize=(12, 6))
        visualize_grid(ax[0], deformed_grid_x1.reshape(N_dim, N_dim), deformed_grid_x2.reshape(N_dim, N_dim))
        visualize_image(ax[1], moving_resampled.reshape(N_dim, N_dim), cmap='gray')
        ax[0].set_title(f"Deformed Grid at Iteration {iteration + 1}")
        ax[1].set_title(f"Resampled Moving Image at Iteration {iteration + 1}")
        plt.tight_layout()
        plt.show()

# Plot energy history
plt.plot(energy_history)
plt.xlabel("Iteration")
plt.ylabel("Energy")
plt.title("Energy vs. Iteration")
plt.show()

In this task, Gauss-Newton optimization was implemented to minimize the sum-of-squared-differences (SSD) between a fixed image and a deformed moving image. Starting with zero weight vectors w1 and w2 (representing no deformation), the algorithm iteratively updated the weights to deform the grid and resample the moving image. At each iteration, the residual error between the fixed and resampled moving image was computed, and the Jacobian matrix was used to estimate how changes in the weights would reduce this error. The weights were then updated to minimize the residual, and the energy (SSD) was tracked across iterations.

The process was repeated for 100 iterations, and at the first iteration and every 20th iteration, the deformed grid and the resampled moving image were visualized. The energy plot should show how the energy decreases over time, indicating that the algorithm was converging toward a solution that better aligned the moving image with the fixed image but in our case it increases and we were not able to figure out why.

Task 5: Experiment with different numbers of Basis Functions

Re-run the registration when the number of basis functions is altered (e.g., $3x3$ or $7x7$ instead of $5x5$) and inspect the outcomes. In your report, comment on differences in the results.

# Running the registration for 3, 5, and 7 basis functions
N_dim = 256  # Image dimension
M_dims = [3, 7]  # Number of basis functions in each direction
delta = 0.001  # Small step for finite differences
iterations = 100  # Number of iterations

for M_dim in M_dims:
    # Parameters
    n_bsplines = M_dim * M_dim   # Number of central basis functions

    # Initialize weight vectors
    w1 = np.zeros(n_bsplines)
    w2 = np.zeros(n_bsplines)

    Phi, phi1, phi2 = generate_splines(N_dim, M_dim)

    # Initialize energy array to store energy at each iteration
    energy_history = []

    # Gauss-Newton optimization loop
    for iteration in range(iterations):
        # Compute the mapped coordinates
        delta_x1 = delta_d(w1, Phi)
        delta_x2 = delta_d(w2, Phi)
        deformed_grid_x1 = grid_x[:, 0] + delta_x1
        deformed_grid_x2 = grid_x[:, 1] + delta_x2
        
        # Resample the moving image using the mapped coordinates
        moving_resampled = map_coordinates(moving, [deformed_grid_x1, deformed_grid_x2], order=1)

        # Compute the residual (difference between fixed and resampled moving image)
        residual = fixed.ravel() - moving_resampled
        
        # Compute the energy (sum of squared differences)
        energy = np.sum(residual**2)
        energy_history.append(energy)

        # Compute Jacobian matrix as described in the slides
        G1 = (map_coordinates(moving, [deformed_grid_x1 + delta, deformed_grid_x2], order=1) - moving_resampled) / delta
        G2 = (map_coordinates(moving, [deformed_grid_x1, deformed_grid_x2 + delta], order=1) - moving_resampled) / delta

        # Jacobian: J = [G1 * Phi, G2 * Phi]
        J = np.hstack([G1[:, np.newaxis] * Phi, G2[:, np.newaxis] * Phi])

        # Solve the linear regression problem: (J^T J)^{-1} J^T r
        JTJ = J.T @ J
        JTr = J.T @ residual
        update = np.linalg.solve(JTJ, JTr)
        
        # Update weights
        w1 -= update[:n_bsplines]
        w2 -= update[n_bsplines:]
        
        # Visualization for the first and every 20th iteration
        if iteration == 0 or (iteration + 1) % 20 == 0:
            fig, ax = plt.subplots(1, 2, figsize=(12, 6))
            visualize_grid(ax[0], deformed_grid_x1.reshape(N_dim, N_dim), deformed_grid_x2.reshape(N_dim, N_dim))
            visualize_image(ax[1], moving_resampled.reshape(N_dim, N_dim), cmap='gray')
            ax[0].set_title(f"Deformed Grid at Iteration {iteration + 1}")
            ax[1].set_title(f"Resampled Moving Image at Iteration {iteration + 1}")
            plt.tight_layout()
            plt.show()

    # Plot energy history
    plt.plot(energy_history)
    plt.xlabel("Iteration")
    plt.ylabel("Energy")
    plt.title("Energy vs. Iteration")
    plt.show()

In this task, we investigated the effect of using different numbers of basis functions (3x3, 5x5, and 7x7) for spline deformation during image registration and compared the outcomes. With fewer basis functions (3x3), the deformation grid was less flexible, resulting in smoother, more global deformations. In contrast, using more basis functions (7x7) allowed for highly localized deformations, leading to better alignment of small image regions but increased the risk of overfitting and instability. The 5x5 case provided a balanced result, allowing both global and local deformations while maintaining stability, with steady energy reduction and good image alignment. Overall, the trade-off between flexibility and stability became apparent as the number of basis functions increased, with more flexible grids potentially leading to overfitting or slower convergence without proper parameter tuning. But as is in the prior task the energy should be decreasing instead of increasing, so we must be doing something wrong, and we are not able to detect the error.

Conclusions

This exercise focused on using nonlinear registration to align MRI brain images by reducing intensity differences with the Gauss-Newton optimization method. Different tasks were completed to understand how various approaches to normalization, smoothing, deformation, and the choice of basis functions affect the results. Using the 99th percentile for normalization worked better than normalizing by the maximum intensity because it reduced the effect of outliers and made the intensity values more comparable between the images.

B-spline basis functions were used to smooth the images while keeping the main structures intact. The weights determined how much each basis function affected different parts of the image, helping to reduce noise effectively. Then changing the weights showed how different values caused localized deformations of the image grid. Applying weights in both directions (W1 and W2) created more complex two-way changes, allowing for more flexible adjustments.

Additionally, the Gauss-Newton method generally updated the weights to reduce differences between images, but the process was not always stable, with the energy level unexpectedly increasing at times. This showed some challenges in ensuring reliable optimization. The number of basis functions used had a big impact on the results. Using fewer functions (3x3) led to smoother, overall changes, while more functions (7x7) allowed for finer adjustments but risked instability. The 5x5 setup offered a good middle ground, balancing flexibility with stability.

Overall, the exercise showed how important it is to choose the right methods for normalization, smoothing, and deformation to achieve accurate image alignment.